using System;
using System.Reflection;
using NBody;
using NBody.Analysis;
using NBody.Cosmology;
using DataFile;
using OptionParser;

[assembly: AssemblyTitle ("NBodyDehnen")]
[assembly: AssemblyVersion ("1.0.*")]
[assembly: AssemblyDescription ("Compute the some stuff from Dehnen (2005).")]
[assembly: AssemblyCopyright ("2005 Joseph D. MacMillan")]

public class DehnenOptions : InputOutputGetOptions
{
    [Option ("-b", "Number of bins")]
    public int NumBins;
    
    [Option ("-s", "Bin spacing")]
    public Bin.Spacing Spacing;
    
    [Option ("-c", "Adjust for centre of mass")]
    public bool CentreOfMass;
    
    [Option ("-r", "Output relative values")]
    public bool Relative;
    
    public DehnenOptions()
    {
        NumBins = 100;
        Spacing = Bin.Spacing.Logarithmic;    
    }
}

public class NBodyDehnen
{
    static double r_200 = 1.0;
    
    public static void Main(string[] args)
    {
        DehnenOptions opts = new DehnenOptions();
        opts.ProcessArgs(args, ".dehnen");
                        
        NBodySystem s = NBodySystem.Read(opts.InFile);
        
        if (opts.CentreOfMass)
            s.AdjustForCentreOfMass(1e-3);
        
        if (opts.Relative)
        {
            Model cosmo = Model.SCDM;
            double rho_c = cosmo.Density(cosmo.Redshift(s.Time));
            r_200 = VirialRadius.CalculateR200(s, rho_c);
        }
        
        Bin den = new DensityBin();
        den.BinSpacing = opts.Spacing;
        den.NumBins = opts.NumBins;
        den.Analyse(s);
        
        VelocityBin disp_r = new VelocityBinDispersion();
        disp_r.VelocityDirection = VelocityBin.Direction.Radial;
        disp_r.BinSpacing = opts.Spacing;
        disp_r.NumBins = opts.NumBins;
        disp_r.Analyse(s);
        
        VelocityBin disp_t = new VelocityBinDispersion();
        disp_t.VelocityDirection = VelocityBin.Direction.Theta;
        disp_t.BinSpacing = opts.Spacing;
        disp_t.NumBins = opts.NumBins;
        disp_t.Analyse(s);
        
        VelocityBin disp_p = new VelocityBinDispersion();
        disp_p.VelocityDirection = VelocityBin.Direction.Phi;
        disp_p.BinSpacing = opts.Spacing;
        disp_p.NumBins = opts.NumBins;
        disp_p.Analyse(s);
        
               
        Table data = new Table(den.Data.Length, 3, "Radius, rho / sigma_r^3, rho / sigma_tot^3");
        for (int i = 0; i < data.Length; i++)
        {
            data[i, 0] = den.Data[i, 0];
            data[i, 1] = den.Data[i, 1] / Math.Pow(disp_r.Data[i, 1], 1.5);
            double tot = 1.0/3.0 * (disp_r.Data[i, 1] + disp_t.Data[i, 1] + disp_p.Data[i, 1]);
            data[i, 2] = den.Data[i, 1] / Math.Pow(tot, 1.5);
        }
        
        data.Print(opts.OutFile, new DataFile.Transform[] { RelRad, Math.Log10, Math.Log10 });
    }
    
    static double RelRad(double r)
    {
        return Math.Log10(r / r_200);
    }
}
